/** ANAYLSIS_21.C
 *  Analyse Macros zum F-Praktikumsversuch 21: Driftkammer
 * 
 *  date:   03-FEB-2012
 *  author: Ch.Weinsheimer
 *  email:  weinshec@students.uni-mainz.de
 */

#define WIRE 16
#define TIME 170

#include <iostream>
#include <string>



/// Define Functions for fitting
double poly4(double* vars, double * pars)
{
  double arg = vars[0]-pars[0];
  return (pars[5] + pars[1]*arg + pars[2]*arg*arg + pars[3]*arg*arg*arg + pars[4]*arg*arg*arg*arg);
}

double gaussian(double* vars, double* pars) {
  return (pars[3]+pars[0]*TMath::Gaus(vars[0],pars[1],pars[2]));
}

double wideningFit(double* vars, double* pars) {
  double arg = vars[0]+pars[0];
  return  (pars[1] + pars[2]*arg*arg + pars[3]*arg*arg*arg);
}



/// Read scintillator count rates for plotting and fitting
void scintillator(char* filename)
{
    // File Opening, Reading and Closing
    FILE* file = fopen(filename, "r");
    Double_t position[15], pmt_left[15], pmt_right[15], coinc[15];
    Double_t dPosition[15], dPmt_left[15], dPmt_right[15], dCoinc[15];
    for(Int_t i=0; i < 15; i++){
      fscanf(file, "%lf %lf %lf %lf", &position[i],&pmt_left[i], &pmt_right[i], &coinc[i]);
      dPosition[i] = 0.5; //position error is constant: 0.5cm
      dPmt_left[i] = TMath::Sqrt(pmt_left[i]);
      dPmt_right[i] = TMath::Sqrt(pmt_right[i]);
      dCoinc[i] = TMath::Sqrt(coinc[i]);
    }
    fclose(file);

    // Create TGraphErrors representation
    TGraphErrors * gr_left = new TGraphErrors(15, position, pmt_left, dPosition, dPmt_left );
    gr_left->SetMarkerColor(kBlue);
    gr_left->SetMarkerStyle(20);
    //    gr_left->SetLineColor(kBlue);
    TGraphErrors * gr_right = new TGraphErrors(15, position, pmt_right, dPosition, dPmt_right );
    gr_right->SetMarkerColor(kRed);
    gr_right->SetMarkerStyle(20);
    //    gr_right->SetLineColor(kRed);
    TGraphErrors * gr_coinc = new TGraphErrors(15, position, coinc, dPosition, dCoinc );
    gr_coinc->SetMarkerColor(kBlack);
    gr_coinc->SetMarkerStyle(20);
    //   gr_coinc->SetLineColor(kBlack);

    // Set Styles
    gROOT->SetStyle("Plain");
    gStyle->SetOptStat(0);

    gr_left->SetTitle("PMT counts");
    gr_left->GetXaxis()->SetTitle("position / cm");
    gr_left->GetYaxis()->SetTitle("counts");

    // Create Canvas
    TCanvas * c1 = new TCanvas("c1", "scintillator()", 1024,768);
    c1->cd(1);
    

    // Draw Graphs
    gr_left->Draw("APE");
    gr_right->Draw("PE Same");
    gr_coinc->Draw("PE Same");

}


/// Read data from file and plot a histogram
TH1F* readHistogram(char* filename, int const N, char* title)
{ 
  // Declaration of fields
  TH1F * h1 = new TH1F(title,title, N, 0, N);
  Int_t value[N];
  string dummy;

  // Open file stream for input
  ifstream in;
  in.open(filename);

  // Read data from file and fill histogramm
  for(int i = 0; i < 3; i++) getline(in, dummy);
  for(int i = 0; i < N; i++) {
    in >> value[i];
    h1->Fill(i,value[i]);
  }
  
  // Close file stream
  in.close();

  // Cosmetics
  h1->GetXaxis()->SetTitle("channel");
  h1->GetYaxis()->SetTitle("counts");

  return h1;
}

/// Read Histogramms and draw them in a Canvas
void drawCanvas(char* wire0, char* wire1, char* time0, char* time1) {

  // Read destination file from stdin
  char outName[16];
  cout << "Output filename: "; 
  cin >> outName;

  // Create and read Histogramms
  TH1F * hWire0 = readHistogram(wire0, WIRE, "wire-0");
  TH1F * hWire1 = readHistogram(wire1, WIRE, "wire-1");
  TH1F * hTime0 = readHistogram(time0, TIME, "time-0");
  TH1F * hTime1 = readHistogram(time1, TIME, "time-1");

  // Create Canvas and Draw
  TCanvas * c1 = new TCanvas("c1","Overview", 1024, 768);
  c1->Divide(2,2);

  c1->cd(1);
  hWire0->Draw();

  c1->cd(2);
  hWire1->Draw();

  c1->cd(3);
  hTime0->Draw();

  c1->cd(4);
  hTime1->Draw();

  // Write maximum values and histogramms to file
  TFile * rFile = new TFile(outName,"RECREATE");
  c1->Write();
  rFile->Close();
  
}


/// Read Histogramms and draw them in a Canvas
void drawCanvasAndFit(char* wire0, char* wire1, char* time0, char* time1, int cut) {

  // Decide Fit Range
  Int_t fit1Max = 170;
  if (cut==1) fit1Max = 100;

  // Read source position from stdin
  Double_t pos, max0, max1, dmax0, dmax1;
  char outName[16];
  cout << "Source position: "; 
  cin >> pos;
  cout << "Output filename: "; 
  cin >> outName;

  // Create and read Histogramms
  TH1F * hWire0 = readHistogram(wire0, WIRE, "wire-0");
  TH1F * hWire1 = readHistogram(wire1, WIRE, "wire-1");
  TH1F * hTime0 = readHistogram(time0, TIME, "time-0");
  TH1F * hTime1 = readHistogram(time1, TIME, "time-1");

  // Create Canvas and Draw
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(111);
  TCanvas * c1 = new TCanvas("c1","Overview", 1024, 768);
  c1->Divide(2,2);

  c1->cd(1);
  hWire0->Draw();

  c1->cd(2);
  hWire1->Draw();

  c1->cd(3);
  TF1 * landau0 = new TF1("landau0", "landau", 0,fit1Max);
  hTime0->Fit("landau0","R");
  max0 = landau0->GetParameter(1);
  dmax0 = landau0->GetParError(1);
  hTime0->Draw();

  c1->cd(4);
  TF1 * landau1 = new TF1("landau1", "landau",0,170);
  hTime1->Fit("landau1","R");
  max1 = landau1->GetParameter(1);
  dmax1 = landau1->GetParError(1);
  hTime1->Draw();
  
  // Write maximum values and histogramms to file
  TFile * rFile = new TFile(outName,"RECREATE");
  c1->Write();
  TNtuple * outTuple = new TNtuple("ntuple", "maximum values from time histogramms", "pos:max0:max1:dmax0:dmax1");
  outTuple->Fill(pos,max0,max1,dmax0,dmax1);
  rFile->Write();
  rFile->Close();

}


/// Read calib-trees in a TChain an distribution
void calibFittedTrees()
{

  // Create TChain and add root files
  TChain chain("ntuple");
  chain.Add("calib_20.root");
  chain.Add("calib_21.root");
  chain.Add("calib_22.root");
  chain.Add("calib_23.root");
  chain.Add("calib_24.root");
  chain.Add("calib_25.root");
  chain.Add("calib_26.root");
  chain.Add("calib_27.root");
  chain.Add("calib_28.root");
  chain.Add("calib_29.root");
  chain.Add("calib_30.root");
  chain.Add("calib_31.root");
  chain.Add("calib_32.root");
  chain.Add("calib_33.root");
  chain.Add("calib_34.root");
  
  chain.Scan();

  Float_t pos[15],max0[15],max1[15],dpos[15],dmax0[15],dmax1[15];
  
  for (int i = 0; i < 15; i++) {
    dpos[i] = 0.05;
    chain.SetBranchAddress("pos", &pos[i]);
    chain.SetBranchAddress("max0", &max0[i]);
    chain.SetBranchAddress("max1", &max1[i]);
    chain.SetBranchAddress("dmax0", &dmax0[i]);
    chain.SetBranchAddress("dmax1", &dmax1[i]);
    chain.GetEntry(i);
  }

  // Create GraphErrors
  TGraphErrors * gY0 = new TGraphErrors(15,pos,max0,dpos,dmax0);
  gY0->GetXaxis()->SetTitle("horizontal position / mm");
  gY0->GetYaxis()->SetTitle("max drifttime");
  gY0->SetTitle("Drift Time vs. Horizontal Position");
  gY0->SetMarkerStyle(20);
  gY0->SetMarkerColor(kRed);

  TGraphErrors * gY1 = new TGraphErrors(15,pos,max1,dpos,dmax1);
  gY1->SetMarkerStyle(20);
  gY1->SetMarkerColor(kBlue);


  // Create Fit Funktion
  TF1 fitY0("fitY0",poly4,2.0,3.4,6);
  fitY0.SetParNames("t0","c1","c2","c3","c4","off");
  fitY0.SetParameters(1.40, 463.7, -580.5, 296.3, -52.6, -70.9);
  fitY0.SetLineColor(kRed);
  gY0->Fit(&fitY0);
  cout << "Chi^2/ndf: " << fitY0.GetChisquare() << " / " << fitY0.GetNDF() << endl;

  TF1 fitY1("fitY1",poly4,2.0,3.4,6);
  fitY1.SetParNames("t0","c1","c2","c3","c4","off");
  fitY1.SetParameters(-2.1, 1164.2, -545.1, 104.0, -7.0, -749.9);
  fitY1.SetLineColor(kBlue);
  gY1->Fit(&fitY1);
  cout << "Chi^2/ndf: " << fitY1.GetChisquare() << " / " << fitY1.GetNDF() << endl;


  TCanvas * c1 = new TCanvas("c1", "c1", 1024,768);
  c1->cd(1);
  gY0->GetYaxis()->SetRangeUser(45.0,65.0);

  gY0->Draw("APE");
  gY1->DrawClone("Same");

}



/// Background calibration estimation
void background(char * wireFilename, char * timeFilename)
{

  // Daclare local fields
  double k, dk;

  // Create new Histogram for evaluation
  TH1F * hEval = new TH1F("hEval", "Integrated Weighted Time Histogram", TIME, 0, TIME);

  // Read histogramms from files
  TH1F * hWire = readHistogram(wireFilename, WIRE, "wire");
  TH1F * hTime = readHistogram(timeFilename, TIME, "time");

  // Calculate average counts from wire histogram
  double sumksqr = 0.0;
  double sumk = 0.0;
  for (int i = 0; i < hWire->GetNbinsX(); i++) {
    double temp = hWire->GetBinContent(i+1);
    sumk += temp;
    sumksqr += temp*temp;
  }
  k = sumk/(16.0*16.0);
  dk = TMath::Sqrt(  (sumksqr-sumk*sumk/16.0)/15.0  )/16.0;
  cout << "average counts: " << k << " +/- " << dk << endl;

  // Integrate time histogramm and multiply with 1/k
  for (int i = 0; i < hTime->GetNbinsX(); i++) {
    hEval->Fill(i,hTime->Integral(0,i)/(k*k));
    double binError = (1/k) * TMath::Sqrt( (dk*hTime->Integral(0,i)/k)*(dk*hTime->Integral(0,i)/k) + hTime->Integral(0,i)  )/(i+1);
    hEval->SetBinError(i+1,binError);
  }

  // Canvas and style
  TCanvas * c1 = new TCanvas("c1", "c1", 1024 , 768);
  gStyle-> SetOptStat(0);
  c1->cd(1);

  hEval->GetXaxis()->SetTitle("drift time / ns");
  hEval->GetYaxis()->SetTitle("distance / cm");
  

  hEval->Draw();

}


/// Analysis of beam widening
void widening(char * file0, char * file1, char * file2, char * file3, char * file4, char * file5, char * file6){

  gStyle->SetOptStat(0);
  gStyle->SetOptFit(111);

  // Create histograms
  TH1F * h0 = readHistogram(file0,WIRE,"wire-82mm");
  TH1F * h1 = readHistogram(file1,WIRE,"wire-85mm");
  TH1F * h2 = readHistogram(file2,WIRE,"wire-88mm");
  TH1F * h3 = readHistogram(file3,WIRE,"wire-91mm");
  TH1F * h4 = readHistogram(file4,WIRE,"wire-94mm");
  TH1F * h5 = readHistogram(file5,WIRE,"wire-97mm");
  TH1F * h6 = readHistogram(file6,WIRE,"wire-100mm");

  // Create Fit Functions
  TF1 f0("f0",gaussian,0.0,16.0,4);
  f0.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f0.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f0.SetNpx(500);
  TF1 f1("f1",gaussian,0.0,16.0,4);
  f1.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f1.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f1.SetNpx(500);
  TF1 f2("f2",gaussian,0.0,16.0,4);
  f2.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f2.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f2.SetNpx(500);
  TF1 f3("f3",gaussian,0.0,16.0,4);
  f3.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f3.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f3.SetNpx(500);
  TF1 f4("f4",gaussian,0.0,16.0,4);
  f4.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f4.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f4.SetNpx(500);
  TF1 f5("f5",gaussian,0.0,16.0,4);
  f5.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f5.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f5.SetNpx(500);
  TF1 f6("f6",gaussian,0.0,16.0,4);
  f6.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f6.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f6.SetNpx(500);
 
  // Do fitting
  h0->Fit(&f0,"R");
  h1->Fit(&f1,"R");
  h2->Fit(&f2,"R");
  h3->Fit(&f3,"R");
  h4->Fit(&f4,"R");
  h5->Fit(&f5,"R");
  h6->Fit(&f6,"R");

  TCanvas * ct = new TCanvas("ct", "ct", 707, 1000);
  ct->Divide(2,4);
  ct->cd(1);
  h0->GetYaxis()->SetRangeUser(0,3000);
  h0->Draw();
  ct->cd(2);
  h1->GetYaxis()->SetRangeUser(0,3000);
  h1->Draw();
  ct->cd(3);
  h2->GetYaxis()->SetRangeUser(0,3000);
  h2->Draw();
  ct->cd(4);
  h3->GetYaxis()->SetRangeUser(0,3000);
  h3->Draw();
  ct->cd(5);
  h4->GetYaxis()->SetRangeUser(0,3000);
  h4->Draw();
  ct->cd(6);
  h5->GetYaxis()->SetRangeUser(0,3000);
  h5->Draw();
  ct->cd(7);
  h6->GetYaxis()->SetRangeUser(0,3000);
  h6->Draw();


  
  // Read Fit Parameters
  double sigma[7], dsigma[7];
  double pos[7], dpos[7];

  sigma[0] = f0.GetParameter(2);
  dsigma[0] = f0.GetParError(2);
  sigma[1] = f1.GetParameter(2);
  dsigma[1] = f1.GetParError(2);
  sigma[2] = f2.GetParameter(2);
  dsigma[2] = f2.GetParError(2);
  sigma[3] = f3.GetParameter(2);
  dsigma[3] = f3.GetParError(2);
  sigma[4] = f4.GetParameter(2);
  dsigma[4] = f4.GetParError(2);
  sigma[5] = f5.GetParameter(2);
  dsigma[5] = f5.GetParError(2);
  sigma[6] = f6.GetParameter(2);
  dsigma[6] = f6.GetParError(2);

  
  // Set pos values
  double val = 8.2;
  for (int i = 0; i < 7; i++) {
    pos[i] = val;
    dpos[i] = 0.05;
    val += 0.3;
    dsigma[i] = 2*sigma[i]*dsigma[i];
    sigma[i] = sigma[i]*sigma[i];
  }
  
  // Create TGraphError
  TGraphErrors * g0 = new TGraphErrors(7,pos,sigma,dpos,dsigma);
  g0->GetXaxis()->SetTitle("d / cm");
  g0->GetYaxis()->SetTitle("#sigma_{gesamt}^{2}");
  g0->SetTitle("Beam Widening Y0");

  // Create Fit function
  TF1 fG0Fit("fG0Fit",wideningFit,8.2,10.0,4);
  fG0Fit.SetParNames("s0","#sigma_{HDC}^{2}","#sigma_{Kollimator}^{2}", "#sigma_{Luft}^{2}");
  fG0Fit.SetParameters(1.1,-5.6,0.16,-0.012);

  g0->Fit(&fG0Fit,"R");

  // Create Canvas
  TCanvas * c0 = new TCanvas("c0","c0",1024,768);
  c0->cd(1);

  g0->Draw("APE");

}



/// Analysis of beam widening
void HVvar1(char * file0, char * file1, char * file2, char * file3, char * file4, char * file5, char * file6){

  gStyle->SetOptStat(0);
  gStyle->SetOptFit(111);

  // Create histograms
  TH1F * h0 = readHistogram(file0,WIRE,"wire - 2.55 kV");
  TH1F * h1 = readHistogram(file1,WIRE,"wire - 2.60 kV");
  TH1F * h2 = readHistogram(file2,WIRE,"wire - 2.65 kV");
  TH1F * h3 = readHistogram(file3,WIRE,"wire - 2.70 kV");
  TH1F * h4 = readHistogram(file4,WIRE,"wire - 2.75 kV");
  TH1F * h5 = readHistogram(file5,WIRE,"wire - 2.80 kV");
  TH1F * h6 = readHistogram(file6,WIRE,"wire - 2.85 kV");

  // Create Fit Functions
  TF1 f0("f0",gaussian,0.0,16.0,4);
  f0.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f0.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f0.SetNpx(500);
  TF1 f1("f1",gaussian,0.0,16.0,4);
  f1.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f1.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f1.SetNpx(500);
  TF1 f2("f2",gaussian,0.0,16.0,4);
  f2.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f2.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f2.SetNpx(500);
  TF1 f3("f3",gaussian,0.0,16.0,4);
  f3.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f3.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f3.SetNpx(500);
  TF1 f4("f4",gaussian,0.0,16.0,4);
  f4.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f4.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f4.SetNpx(500);
  TF1 f5("f5",gaussian,0.0,16.0,4);
  f5.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f5.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f5.SetNpx(500);
  TF1 f6("f6",gaussian,0.0,16.0,4);
  f6.SetParNames("Constant", "Mean", "Sigma", "Offset");
  f6.SetParameters(2000.0, 8.5, 0.5, 50.0);
  f6.SetNpx(500);
 
  // Do fitting
  h0->Fit(&f0,"R");
  h1->Fit(&f1,"R");
  h2->Fit(&f2,"R");
  h3->Fit(&f3,"R");
  h4->Fit(&f4,"R");
  h5->Fit(&f5,"R");
  h6->Fit(&f6,"R");

  TCanvas * ct = new TCanvas("ct", "ct", 707, 1000);
  ct->Divide(2,4);
  ct->cd(1);
  h0->GetYaxis()->SetRangeUser(0,3000);
  h0->Draw();
  ct->cd(2);
  h1->GetYaxis()->SetRangeUser(0,2000);
  h1->Draw();
  ct->cd(3);
  h2->GetYaxis()->SetRangeUser(0,5000);
  h2->Draw();
  ct->cd(4);
  h3->GetYaxis()->SetRangeUser(0,2000);
  h3->Draw();
  ct->cd(5);
  h4->GetYaxis()->SetRangeUser(0,2000);
  h4->Draw();
  ct->cd(6);
  h5->GetYaxis()->SetRangeUser(0,3000);
  h5->Draw();
  ct->cd(7);
  h6->GetYaxis()->SetRangeUser(0,2000);
  h6->Draw();


  
  // Read Fit Parameters
  double sigma[7], dsigma[7];
  double pos[7], dpos[7];

  sigma[0] = f0.GetParameter(0);
  dsigma[0] = f0.GetParError(0);
  sigma[1] = f1.GetParameter(0);
  dsigma[1] = f1.GetParError(0);
  sigma[2] = f2.GetParameter(0);
  dsigma[2] = f2.GetParError(0);
  sigma[3] = f3.GetParameter(0);
  dsigma[3] = f3.GetParError(0);
  sigma[4] = f4.GetParameter(0);
  dsigma[4] = f4.GetParError(0);
  sigma[5] = f5.GetParameter(0);
  dsigma[5] = f5.GetParError(0);
  sigma[6] = f6.GetParameter(0);
  dsigma[6] = f6.GetParError(0);

  
  // Set pos values
  double val = 2.55;
  for (int i = 0; i < 7; i++) {
    pos[i] = val;
    dpos[i] = 0.01;
    val += 0.05;
    //    dsigma[i] = 2*sigma[i]*dsigma[i];
    //    sigma[i] = sigma[i]*sigma[i];
  }
  
  // Create TGraphError
  TGraphErrors * g0 = new TGraphErrors(7,pos,sigma,dpos,dsigma);
  g0->GetXaxis()->SetTitle("cathode voltage / kV");
  g0->GetYaxis()->SetTitle("Constant");
  g0->SetTitle("Variation Of Cathode Voltage Constant Fraction");

  /*
  // Create Fit function
  TF1 fG0Fit("fG0Fit",wideningFit,8.2,10.0,4);
  fG0Fit.SetParNames("s0","#sigma_{HDC}^{2}","#sigma_{Kollimator}^{2}", "#sigma_{Luft}^{2}");
  fG0Fit.SetParameters(1.1,-5.6,0.16,-0.012);

  g0->Fit(&fG0Fit,"R");
  */

  // Create Canvas
  TCanvas * c0 = new TCanvas("c0","c0",1024,768);
  c0->cd(1);

  g0->Draw("APE");

}
